{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using JuLIP, JuLIP.ASE, JuLIP.Potentials \n",
    "using PyPlot\n",
    "\n",
    "function plot2d(at, ttl)\n",
    "    X = positions(at) |> mat \n",
    "    plot(X[1,:], X[2,:], \"b.\", markersize=20)\n",
    "    title(ttl)\n",
    "    PyPlot.draw()\n",
    "    PyPlot.pause(0.0001)    \n",
    "end \n",
    "\n",
    "# For now, this notebook is just to implement a temporary\n",
    "# example for MD, Newtonian dynamics. This functionality\n",
    "# should be moved into the JuLIP library, using more\n",
    "# sophisticated integrators. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "at = bulk(\"Si\", cubic=true) * (4,4,1)\n",
    "set_pbc!(at, true)\n",
    "set_calculator!(at, StillingerWeber())\n",
    "set_constraint!(at, FixedCell(at));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this should be done properly, using in particular \n",
    "# a proper scaling of the gaussian depending on temperature \n",
    "q = position_dofs(at) \n",
    "p = randn(length(q)) * 0.1\n",
    "set_momentum_dofs!(at, p) \n",
    "write(\"md.xyz\", at)\n",
    "# prepare for plotting\n",
    "pygui(true)\n",
    "PyPlot.ion()\n",
    "PyPlot.hold(false)\n",
    "# create a counter to save only every 10 steps \n",
    "cnt = 0 \n",
    "dt = 0.003\n",
    "E0 = energy(at, q)\n",
    "Ep = [energy(at, q)-E0]\n",
    "Ek = [0.5 * norm(p)^2] \n",
    "for n = 1:10_000\n",
    "    p -= dt * gradient(at, q)    # Euler-A method\n",
    "    q += dt * p\n",
    "    cnt += 1\n",
    "    push!(Ep, energy(at, q) - E0)\n",
    "    push!(Ek, 0.5 * vecnorm(p)^2)\n",
    "    if cnt == 40 \n",
    "        cnt = 0 \n",
    "        # write current configuration into at\n",
    "        set_position_dofs!(at, q) \n",
    "        set_momentum_dofs!(at, p)\n",
    "        # write to file\n",
    "        write(\"md.xyz\", at, :append)\n",
    "        plot2d(at, \"n = $(n), Ep = $(Ep[end]), Ek = $(Ek[end])\") \n",
    "    end \n",
    "end \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PyPlot.pygui(false)\n",
    "x = collect(1:length(Ep)) * dt\n",
    "plot(x, Ep, \"b-\", x, Ek, \"g-\", x, Ep+Ek, \"r-\")\n",
    "legend((\"E_pot\", \"E_kin\", \"E_tot\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
